import numpy as np
from scipy.optimize import fsolve
from scipy import integrate
from scipy.stats import norm
import seaborn as sns
import matplotlib.pyplot as plt
from tqdm import tqdm
import math


# 以beta分布的参数alpha和beta为未知数的方程组
# 该方程组是根据高度均值和变异系数列出的
# 左跨方程组
def func_left(param_list):
    # x, y 即 alpha, beta
    x, y = param_list[0], param_list[1]
    return [x / (x + y) - h_l_mean,
            x * y / ((x + y) ** 2 * (x + y + 1) * h_l_mean ** 2) - Cv_l ** 2]


# 右跨方程组
def func_right(param_list):
    # x, y 即 alpha, beta
    x, y = param_list[0], param_list[1]
    return [x / (x + y) - h_r_mean,
            x * y / ((x + y) ** 2 * (x + y + 1) * h_r_mean ** 2) - Cv_r ** 2]


# 生成两跨高度的随机数, 各num个
def generate_height(num):
    ans_left = fsolve(func_left, np.array([1, 1]))
    ans_right = fsolve(func_right, np.array([1, 1]))
    alpha_left, beta_left = ans_left[0], ans_left[1]
    alpha_right, beta_right = ans_right[0], ans_right[1]
    h_l_rand = np.random.beta(a=alpha_left, b=beta_left, size=num)
    h_r_rand = np.random.beta(a=alpha_right, b=beta_right, size=num)
    return h_l_rand, h_r_rand


def solve_X1(I_l, I_r, P):
    """
    # 1处的支座反力X_1. 公式推导见word文档2.1.1
    :param I_l: 左跨惯性矩
    :param I_r: 右跨惯性矩
    :param P: P可为外荷载或单位荷载
    :return: 由C点外荷载或单位荷载引起的X1
    """
    Delta_1P = -1 / (E * I_l) * P * L ** 3 / 6 - 1 / (E * I_r) * 5 * P * L ** 3 / 48
    delta_11 = L ** 3 / 3 * (1 / E / I_l + 1 / E / I_r)
    X1 = -Delta_1P / delta_11
    return X1


def M(x, I_l, I_r, P):
    """
    计算M_bar和M_P. 计算公式见word文档2.1.2
    :param L: 每跨长
    :param x: 位置参数
    :param I_l: 左跨惯性矩
    :param I_r: 右跨惯性矩
    :param P: 外荷载或单位荷载
    :return: M_bar or M_P
    """
    X_1 = solve_X1(I_l, I_r, P=P)
    F_B = 3 * P / 2 - 2 * X_1
    F_A = X_1 - P / 2
    if x < 0:
        print("x out of range")
        return
    elif x <= L:
        return F_A * x
    elif x <= 3 * L / 2:
        return F_A * x + F_B * (x - L)
    elif x <= 2 * L:
        return F_A * x + F_B * (x - L) - P * (x - 3 * L / 2)
    else:
        print("x out of range")
        return


# 两个弯矩以及抗弯刚度的乘积, 分0-L、L-3L/2、3L/2-2L三段
# 0-L
def Mp_Mbar_1(x, I_l, I_r):
    return M(x, I_l, I_r, P_C) * M(x, I_l, I_r, 1) / (E * I_l)


# L-3L/2
def Mp_Mbar_2(x, I_l, I_r):
    return M(x, I_l, I_r, P_C) * M(x, I_l, I_r, 1) / (E * I_r)


# 3L/2-2L
def Mp_Mbar_3(x, I_l, I_r):
    return M(x, I_l, I_r, P_C) * M(x, I_l, I_r, 1) / (E * I_r)


def plot_pdf(data_list):
    """
    绘制概率密度函数
    :param data_list: Delta_C的取值列表
    :return: None
    """
    sns.displot(data_list, kind="kde", bw_adjust=0.25)
    plt.xlabel("Delta_C")
    plt.show()


def calculate():
    """
    进行所有所需结果的计算
    :return: a tuple: (失效概率的估计值, C处位移值列表)
    """
    fail_num = 0
    h_l_random, h_r_random = generate_height(n)  # 生成高度列表
    print("高度随机值生成完毕!")
    I_l_list = [breadth * h ** 3 / 12 for h in h_l_random]  # 左跨惯性矩的列表
    I_r_list = [breadth * h ** 3 / 12 for h in h_r_random]  # 右跨惯性矩的列表
    print("\n惯性矩随机值生成完毕!\n下面开始进行{}次模拟!".format(n))
    Delta_C_ls = []
    for i in tqdm(range(n)):
        I_l = I_l_list[i]
        I_r = I_r_list[i]
        tuple_1 = integrate.quad(Mp_Mbar_1, 0, L, args=(I_l, I_r))
        tuple_2 = integrate.quad(Mp_Mbar_2, L, 3 * L / 2, args=(I_l, I_r))
        tuple_3 = integrate.quad(Mp_Mbar_3, 3 * L / 2, 2 * L, args=(I_l, I_r))
        Delta_C = tuple_1[0] + tuple_2[0] + tuple_3[0]
        if Delta_C > Delta_C_max:
            fail_num += 1
        Delta_C_ls.append(Delta_C)
    Pf_hat = fail_num / n  # 失效概率估计值
    Var_Pf_hat = Pf_hat * (1 - Pf_hat) / n  # 失效概率估计值的方差估计值
    V_Pf_hat = math.sqrt((1 - Pf_hat) / (n * Pf_hat))  # 失效概率估计值的变异系数
    u = norm.ppf(1 - alpha / 2)  # 标准正态分布的上α/2分位数
    CI_left = Pf_hat - u * math.sqrt(Var_Pf_hat)  # 置信区间左侧边界
    CI_right = Pf_hat + u * math.sqrt(Var_Pf_hat)  # 置信区间右侧边界
    epsilon_max = u * V_Pf_hat  # 失效概率估计值的相对误差的上限
    print("{}次模拟下失效概率为: {}".format(n, Pf_hat))
    print("总体参数失效概率的置信区间为[{}, {}]".format(CI_left, CI_right))
    print("失效概率估计值的相对误差的上限为: {}".format(epsilon_max))
    return Pf_hat, Delta_C_ls


if __name__ == '__main__':
    # 设置初始参数
    np.random.seed(2021)  # 设置随机数种子
    n = 200000  # 模拟次数
    P_C = 50 * 1e3  # C处外荷载(单位为N)
    L = 6  # 每跨跨度(m)
    breadth = 0.2  # 跨的宽度(m)
    E = 28000 * 1e6  # 弹性模量(Pa)
    h_l_mean = 0.45  # 左胯高度均值(m)
    Cv_l = 0.1  # 左胯高度变异系数
    h_r_mean = 0.4  # 右胯高度均值(m)
    Cv_r = 0.05  # 右胯高度变异系数
    Delta_C_max = 0.0042  # 右胯跨中位移限值(m)
    alpha = 0.05  # 显著性水平

    # 原问题计算
    data = calculate()[1]

    # 画位移的概率密度函数图
    plot_pdf(data)

    # 改变外荷载值, 计算对应失效概率
    print("************************************************")
    print("以下开始计算外荷载改变的情况")
    n = 10000  # 设置每个荷载情况下的模拟次数
    P_C_ls = []  # 横坐标: 外荷载
    Pf_ls = []  # 纵坐标: 失效概率
    i = 0
    for P_C in np.arange(30 * 1e3, 55 * 1e3, 1e3):
        Pf_ls.append(calculate()[0])
        P_C_ls.append(P_C)
        i += 1
        print("******************************")
        print("第{}个外荷载值情况下的计算完毕".format(i))
        print("******************************")
    sns.set_theme(style="darkgrid")
    sns.lineplot(x=P_C_ls, y=Pf_ls)
    plt.ylabel("Pf")
    plt.xlabel("Pc(N)")
    plt.show()
